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The spectrum of two bosons and two fermions in a trap is calculated using a correlated-Gaussian 
basis throughout the range of a broad Fano-Feshbach resonance. The calculations provide a few- 
body solution to the magneto-association of fermionic Feshbach molecules. This solution is used 
to study the time evolution of the system as the scattering length changes, mimicking experiments 
with Bose-Fermi mixtures near Fano-Feshbach resonances. The structure of avoided crossings in 
the few-body spectrum enables an interpretation of the dynamics of the system as a sequence of 
Landau-Zener transitions. The calculated molecule formation rate is compared with experimental 
observations. 



I. INTRODUCTION 

Cold atomic gases are of interest to physicists because 
the fundamental behavior of quantum matter can be in- 
vestigated in such systems. After Bose-Einstein conden- 
sates p] and degenerate Fermi gases [2] were created, 
Fano-Feshbach resonances were used to tune the inter- 
action between atoms and to make molecules at ultra- 
cold temperatures. Since those early developments, these 
resonances have become essential for manipulating ul- 
tracold gases [21 H] ■ Fano-Feshbach resonances in Bose- 
Fermi mixtures sparked research interest later [5J [3] , be- 
cause of the rich internal structure and the anisotropic 
interaction of heteronuclear Feshbach molecules. Load- 
ing fermionic heteronuclear molecules into optical lat- 
tices reduces their inelastic collisions and extends their 
lifetimes [7lU2j. These experiments have motivated the 
development of theories to describe the dynamical pro- 
cesses of ultracold gases that occur during field ramps 
across a Fano-Feshbach resonance, which have been car- 
ried out from the perspectives of both many-body and 
few-body physics, for Bose-Einstein condensates and two- 
component Fermi gases [T3J 114] . However it has been 
challenging to describe the magneto-association of Fesh- 
bach molecules formed from a Bose-Fermi mixture. 

This article explores the spectrum and dynamics of 
three or four trapped atoms and shows how a few-body 
formulation can obtain accurate solutions of the system 
without making the standard approximations of many- 
body theory. This formulation provides an explicit rep- 
resentation of the avoided crossings responsible for inter- 
conversion between atomic states in a Bose-Fermi mix- 
ture and fermionic molecular states. Our results directly 
apply to optical lattice experiments with a few parti- 
cles on one site, and could also be applied to larger 
Bose-Fermi mixtures through a frequency rescaling as in 
[P31 ITS] , The calculated dynamics provide a few-body 
perspective on the process of magneto-association of Fes- 
hbach molecules as varying degrees of adiabaticity and 
diabaticity. 

Specifically, to describe the behavior of an ultracold 
Bose-Fermi mixture from a few-body perspective in the 
vicinity of a Fano-Feshbach resonance, the spectrum is 
computed for two systems: one consisting of two bosons 



and one fermion (BBF) and the other consisting of two 
bosons and two fermions (BBFF). Although the one bo- 
son and two fermion system is also a subsystem of the 
four body system(BBFF), since it will not provide in- 
formation of the trimer, we did not explore its property. 
In each system, interactions occur only through a short- 
range potential. To carry out this study for a concrete ex- 
ample, we take 87 Rb as the boson and 40 K as the fermion. 
The s-wave scattering length a characterizing the short- 
range interaction between a boson and a fermion pair 
is tuned in the standard manner |16| . This model al- 
lows us to explore the magneto-association of fermionic 
molecules near a broad Fano-Feshbach resonance. Our 
solution of the few-body problem in this universal regime 
gives various properties accurately, including energy lev- 
els and the ramped time-dependent dynamics of the full 
quantum mechanical system. This type of analysis en- 
ables a detailed understanding of the global topology of 
the spectrum. The transition probabilities of this system 
into alternative possible final configurations are quantita- 
tively predicted, for the case in which interactions change 
with time, as in field-ramp experiments |12j . 



II. METHOD 

Eigenfunctions and eigenenergies for the few-body 
Hamiltonian [Eq.([T)] are determined using a variational 
calculation that utilizes a correlated-Gaussian basis [T7] . 
This method has been previously used to treat the two- 
component fermion system |13j and identical boson sys- 
tem |18| . After determination of the adiabatic spectrum, 
the time-dependent Schrodinger equation is solved nu- 
merically using the diabatic-by-sector method [331 [T9] , 
The few-body Hamiltonian adopted has only pairwise in- 
teractions between bosons and fermions in an isotropic 
harmonic trap: 
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where B denotes boson, F denotes fermion, and the inter- 
particle potential is V = Vq exp(— ), i.e., an attractive 
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Gaussian. The width do of the Gaussian is fixed, and 
the depth Vq is tuned to produce the desired two-body 
scattering length a. The boson-fermion reduced mass /i, 



defined as fi = 



- , is half of the mass unit in our cal- 
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culation. To obtain results independent of the details of 
the model potential, the range of the potential is limited 
to c?o "C a-ho, where = (h/2nu} Q ) 1 / 2 is the trap length 
based on this reduced mass. In the present calculation, 
do = O.Ola/jo is used. The forms of the wave functions 
represented by this correlated Gaussian basis are listed 
in Table U 

Each basis functions defined in Table U is character- 
ized by a set of average interparticle distances and 
by the orientations of Jacobi vectors itj (the Ui are not 
needed for total angular momentum L = case). These 
entities are selected semi-randomly from several typical 
and intermediate configuration functions. Since both the 
long-range and short-range portion of the wave functions 
require an accurate description, the dij range from a frac- 
tion of do up to several a^ . The typical size of the total 
basis set used in these calculations varies from 1000 to 
5000. The advantage of the correlated- Gaussian basis set 
is that all the matrix elements can be evaluated analyti- 
cally. One disadvantage of the correlated- Gaussian basis 
is that a large basis set may have a strong linear depen- 
dence, possibly resulting in numerical instability. Thus, 
we apply a method to suppress linear dependence during 
the generation and optimization process. The basis set 
is then fixed while Vo is tuned to give different two-body 
scattering lengths. Matrix elements are calculated once 
and then used to obtain the spectrum throughout the 
entire range of the Fano-Feshbach resonance. 

The basis functions are tested for their convergence 
and accuracy in describing the system. The states in 
a ~ + and a ~ 0~ limits are in agreement with the 
limiting-analytical behavior. The calculation elucidates 
the spectrum as a function of the scattering length and 
of the ramping dynamics. The inverse scattering length, 
which is the adiabatic parameter A of the ramping dy- 
namics, is calculated in free space for the Gaussian po- 
tential by solving the equation of relative motion between 

two particles : -<£"(r) + exp(-^))0(r) = 0. The 

zero-energy scattering solution generally has an asymp- 
totic form equal to <p(r) —> C(r—a), where C is a constant 
and a is the s-wave scattering length. 

With an optimized basis set, the spectra are calcu- 
lated as a function of l/o, which presents a series of 
avoided crossings between the adiabatic levels near the 
point of trimer's formation and near unitarity(l/a ~ 0). 
Each avoided crossing can be qualitatively characterized 
by its width 8X, and by the range over which at least 
two adiabatic energy levels interact strongly, and grouped 
into two main categories: narrow crossings, whose 5X > 
1/a/ioj an d wide crossings, whose SX » l/a/, . Narrow 
crossings are approximated as exact, uncoupled crossings 
in the diabatization procedure, which gives a smooth, 
physical, partially diabatized spectrum as a function of 
X/a. Only the wide crossings are physically relevant for 



the range of ramping speeds of current interest, so the 
diabatization procedure leaves the wide crossings adia- 
batic. 

Thus the idea of diabatization is to compare the degree 
of similarity (quantitatively, the overlap) of the eigen- 
functions at two nearby interaction strengths, by con- 
trolling the distance between the two points being com- 
pared, the best partially diabatized spectrum is selected. 
This selection is somewhat subjective, but after an appro- 
priate diabatization, the structure of the avoided cross- 
ings permits a global qualitative view of the evolution of 
the system through different pathways from the weakly- 
interacting Bose-Fermi mixture at a — > 0~ to a strongly- 
interacting gas at a — > + . These pathways cross the 
universal region where a goes from — oo to oo, mimick- 
ing the experiments carried out in different laboratories 
[THI2U]. Fig [I] presents the partially diabatic spectrum 
in the resonance region. 

A more quantitative way to describe the avoided cross- 
ings is to evaluate the first derivative P-matrix coupling 
between the two energy eigenstates. The P-matrix is 
the nonadiabatic coupling between two adiabatic states, 



specifically: Pij — ^il^f ), where A is the adiabatic pa- 
rameter. A narrow avoided crossing corresponds to a nar- 
row and sharply peaked P-matrix centered at the cross- 
ing point, while a wide avoided crossings corresponds to 
a wide and smooth P-matrix. The diabatization proce- 
dure could also be carried out by calculating the P-matrix 
elements connecting each pair of adiabatic states and se- 
lecting those that are sufficiently wide according to some 
quantitative criterion to determine the adiabatic states. 

After calculating the spectrum, the initial configura- 
tion is propagated using the time-dependent Hamilto- 
nian. Starting from the ground state of the noninter- 
acting limit (a — > 0~), the parameter A (A = — ) is 
ramped through the resonance to the strongly interact- 
ing limit (a — > + ) at different speeds v = 4|. The 
Landau- Zener approximation is applied to interpret our 
results. The Landau-Zener approximation, predicts that 
the probability for a transition from the adiabatic tyj to 
tyj is Tij = exp(— Xij/x)i where x is the sweep rate of 
the adiabatic variable, and Xij is a parameter extracted 
from properties of the adiabatic eigenstates. 

The nonadiabatic coupling (P-matrix) controls the 
probability of nonadiabatic transitions. Clark has shown 
that if a transitions has a form consistent with Landau- 
Zener approximation, then the P-matrix element for a 
transition from v&j to \tj has a Lorentzian form whose 
width(or height), along with the corresponding eigenen- 
ergies, characterize the Landau-Zener parameter Xij |21j . 
Consequently, we evaluate all the potentially important 
Pij numerically and verify that the couplings between 
low-lying configurations generally have a smooth single- 
peak form that is approximately Lorentzian. However, 
the couplings between high-lying configurations do not 
have a perfect Lorentzian form, and can be multipeaked, 
as shown in Fig. [3] and Table III The most important 
for this specific dynamical sweep corresponds to transi- 
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TABLE I: Possible configurations and Gaussian basis functions. BBF denotes the trimer, BF denotes the dimer, and B (or F) 
denotes a single boson (fermion). R cm is the center-of-mass coordinate, ipo is the ground state of the center-of-mass motion 
in a harmonic trap ^o(Rcm) = e _2iicm ' /a ' i ° , and S denotes all the symmetrization or antisymmetrization appropriate for the 
identical bosonic or fermionic particles, S — (1 + Pb,,bj){l — Pfcjj), and P is the permutation operator. 

System configurations basis function 

BBFL = 0+ BBF, BF+B, B+B+F ^(n, r 2 , r 3 ) = S ^Po^cmK ^>' r « /2d « } 

BBFF L = 0+ BBF+F, BF+B+F, B+B+F+F *d, u (ri, r 2 , r 3 , r 4 ) = S U {TL cm )e- ^» r « /M « '} 

BBFF L = 1- BBF+F, BF+BF, BF+B+F, B+B+F+F * d>u (ri, r 2 , r 3 , r 4 ) = S {^(RcmXEi u * z *) e ~ Ej>1 } 



tions among several configurations are listed in Table [IT] 
The final probability distribution can be explained as a 
sequence of Landau-Zener transitions between these par- 
tially diabatic classes where the positions of the P-matrix 
elements' peaks determine a specific order in which the 
transitions occur. 

Specifically for our system, the formation of the trimer 
happens at a negative scattering length, whereas the 
formation of the dimer(s) happens at unitarity. The 
Landau-Zener parameters obtained from the P-matrix 
analysis are listed in Table |LU} The Landau-Zener model 
shows qualitative agreement with the numerical results. 
Fig. [2] presents the probability distribution of three or 
four particles system with total angular momentum 
and 1 following a unidirectional ramp, as a function of 
the ramping speed 4^. The numerical ramps are initi- 
ated at Xi ~ — 8/a/jo and finalized at A/ ~ 8/a; lo , where 
do <C a sc <C aho is satisfied, i.e. within the universal 
regime. As the speed is increased, the most probable fi- 
nal state changes from a bound trimer configuration to a 
bound dimer configuration and then to an atomic ground 
state, which is consistent with experimental observations. 

To connect our few body calculation to many body ex- 
periments, we construct a dimensionless ramping speed 
which has the density information of the system built in: 
X = yf \ I , in which M is the total mass of the molecule, 
p is the density of boson-fermion pair in the noninter- 
acting system, and A = rf( - 1 ^° ac ' 1 is the ramping speed 
at which the inverse scattering length is ramped. Our 
calculation provides probability distribution into various 
final states across a Feshbach resonance as a function 
of the dimensionless ramping speed. At the same time, 
the dimensionless ramping speed could be extracted from 
experimental data at each value of our calculation 
predicts the molecule formation at that ramping speed. 
The above definition of dimensionless ramping speed is 
in a homogeneous system, however, we propose to use 
the peak density in an inhomogeneous system such as 
ultracold gases in harmonic trap. 

In other words, the connection between the few body 
and many body system is realized by solving a few body 
problem in trap tight enough so that the peak density 
of the few body system is the same as the peak density 
of the many body system. The rescaled oscillator length 



depends on the peak density of the many body system 
a ho — (p/ Pexp) 1 ^ 3 , where p is the dimensionless peak 
density p — p * (a^ Q ) 3 from the few-body calculation, 
Pexp is the peak density measured in experiments. For 
our calculation, the ratio of artificial trapping frequency 
to the experimental trapping frequency is around 
30. 

In our calculation, the following units are used: H = 1 , 
boson-fermion reduced mass a — mB " F = 0.5, the 

mass ratio = fX, total mass of the molecule M = 

rap 4U ' 

rriB + mp = 2.31, the trap frequency u> = 1, the os- 
cillator length of the boson-fermion relative motion is 
= (j^j) 1 ^ 2 = 1' which is set to be the length unit, 
and the time unit — . Density of molecules is defined 
as: J pd 3 r = N. For experimental system, N is around 
10 4 ~ 10 5 [Hll], while for our calculation, N is 2. Di- 
mensionless peak density p — 1.066/a^ o . 

The artificially tight trap could be related to the num- 
ber ratio of the few body system to the many body sys- 
tem. For experimental conditions, the temperature of the 
fermionic atoms is around 0.2 ~ 0.3Tp, and the number 
of particles is around 10 4 to 10 5 [12122], so the density 
profile of the system could be described by the zero tem- 
perature Thomas- Fermi approximation |23j . The peak 
density of a zero temperature noninteracting Fermi gas 
in an isotropic harmonic trap is proportional to N 1 / 2 , 
when the particle number of N is much large than 1. So 
the relation between the experimental trapping frequency 
and the rescaled trapping frequency could be expressed 
in terms of the particle number ratio of few body calcula- 
tion to many body experiments: the ratio ~ (^)^- 
The temperature of bosons in the Bose-Fermi mixture 
is around 1.1 ~ L2T C , thus the boson's peak density is 
tuned to match the peak density of the fermions in ex- 
periments. So the fermion peak density itself could give 
a good estimation of the maximum number of atom pairs 
near the center of the trap. 



III. RESULTS 

To relate our results to recent JILA experiments |12) . 
the Landau-Zener parameter for trimer-dimer-atom tran- 
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sitions 6 — Xmol/x is cas t m terms of experimentally ac- 
cessible variables. Since the trimer and the dimer are 
not distinguished experimentally, we present one analy- 
sis for which the trimer formation is counted and another 
analysis where it is excluded from the count of molecules 
formed. 

The adiabatic parameter in our calculation can be 
related to the ramping experiments, through an as- 
sumption that the dependence of the two-body scat- 
tering length a (B) on the magnetic field is approxi- 
mated in the usual manner as a(B) = a og (l + b-b )' 
whereby = (dB / dt) / (wa^g) . Therefore, S = 

Xmoi(dB / dt)~ x ph\wabg\/m. The resonance is assumed 
to be sufficiently broad that, in the ramping range of 
magnetic field, to B — Bo, whereby the ramping speed 
d ^ 1 l°- s ^ can be simplified to — — 4§ ■ This method allows 

at ^ a,bg uj at 

us to interpret experimental data using /3(4r). 

The molecule fraction in the strongly interacting limit 
is the most relevant quantity to compare with experi- 
ments. For the present 4-body system, the molecular 
fraction is defined as the probability of ending up in 
the dimer-dimer configuration (this only occurs in the 
L = case) plus half of the probability of the BF+B+F 
configuration. If we count each trimer as a molecular 
bound state also, the trimer probability should be added. 



A. Energy spectrum 

The energy spectrum is the first thing to calculate and 
it offers the launching point for further analysis. The sys- 
tems selected for study are motivated by consideration of 
the relevant energy to the experimental temperature un- 
der the angular momentum and parity restriction. Since 
the temperature is in the ultracold regime, only low par- 
tial waves scattering can contribute. The fermionic na- 
ture of the dimer and trimer molecules formed implies 
that the ground state of the non-interacting system and 
the ground state of the strongly interacting system have 
different parities, so both parties are considered. 

The BBF L = + system is a subsystem of the four 
body system (BBFF), and it is also computed because 
it provides the properties of the trimer. The BBFF 
L = + ground state is the ground state of the four- 
particle BBFF system in the a ~ 0+ limit. The BBFF 
L = 1~ ground state is the ground state of the four- 
particle BBFF system in the a ~ CP limit. In this system 
a two-dimer configuration is allowed, thus it is possible to 
enhance the formation rate of the fermionic dimers com- 
pared to the BBFF L = + state. Because the dimers 
formed are identical fermions, the relative angular mo- 
mentum between two s-wave dimers must be odd, which 
makes states with nonzero total angular momentum state 
possibly important. This is different from the bosonic 
dimer case [51 [53] 

In each calculation, the various classes of final states 
are distinguished by their configuration in the a ~ + 



limit. Both the energy and the pair-correlation function 
are examined to determine the state configuration. The 
sublevels within each class have different curvatures near 
avoided crossings. The "sublevels" that are close in en- 
ergy just before reaching the crossing region are often 
degenerate states in the noninteracting limit. However 
the boson-fermion pair in these states have differences in 
their relative motion, e.g., different angular momentum 
or a different number of radial nodes; consequently the 
ramp affects these states differently, especially near the 
avoided crossings. In the a ~ + regime, levels that may 
not be energetically close in the noninteracting limit, can 
sometimes end up in the same configuration, analogous 
to the promotion of diatomic electronic orbitals in the 
united atom versus the separated atom limits. All of the 
present calculations are performed at w = 1, a/ lo = 1, 
m B = 1.5815, m F = 0.729 88 ( 40 K and 87 Rb). Figure [l] 
presents the partially-diabatized energy spectrum for the 
few-body systems considered here. 



B. Time evolution 

The time-evolution of the initial state correlates with 
the adiabatic energy spectrum, e.g., when the relative 
probabilities change rapidly at some time to, there is 
usually a corresponding avoided crossing in the energy 
spectrum at to- But in the weak coupling regime, the 
probability in each class remains steady and the phase of 
the probability amplitude is described by the WKB-type 
evolution. 

The final probability distribution predicted by the nu- 
merical time evolution is fitted to the Landau-Zener 
frame work by including the minimum number of nec- 
essary transitions. The two main criteria for ranking the 
importance of any given transition are (1) the initial pop- 
ulation of the relevant states, indicating the population 
that can be transferred, (i.e., if the initial population of 
two states are both small, even if the coupling between 
these two levels are strong, this transition's contribution 
to the final probability distribution is still minor) and (2) 
the Landau-Zener parameter given by our computed P- 
matrix analysis, characterizing the strength of the transi- 
tion and its corresponding characteristic ramp rate. The 
optimal fitting result is shown in Fig. [2] and Table [IT] 
In general, as the energy spectrum gets more compli- 
cated, the minimum number of transitions needed to get 
a qualitatively good fit to the numerical time-evolution 
increases. 

The fit shows that our numerical results appear to 
be well described by an incoherent Landau-Zener model. 
However, in principle, to analyze the interference of tran- 
sition amplitude pathways leading to indistinguishable fi- 
nal states, the phase information must be included, as in 
the Stueckelberg case [25], where inclusion of the phase 
information was essential to determine the oscillation of 
probabilities as functions of energy. In our calculation, 
the incoherent Landau-Zener model fits adequately, as 
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FIG. 1: The energy spectra are plotted as functions of the dimensionless inverse scattering length in trap oscillator units aho/dsc 
for the three few-body systems considered. Classes and typical sub-levels in each class are pointed by arrows on the figure, (a) 
BBF L — + . (b) BBFF L — + . (c) BBFF L = 1~. The inset shows avoided crossings and a pictorial demonstration of our 
diabatization near unitarity. 



was the case for the two-component fermion system |13j . 
One possible reason is that the initial population for each 
important transition has always been chosen to reside in 
a single state, which results in little if any interference. 



Some evidence exists to support this suspicion; specifi- 
cally, when the initial state chosen for the time evolution 
is picked to be the ground state, a probability oscillation 
is observed in some cases as a function of the ramping 
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FIG. 2: The final probability distribution are plotted versus the ramping speed. Starting from the ground state of the nearly 
noninteracting fewbody system, as the dimensionless ramping speed x — H^l4tl (A = 1 ' ** dt I) ' s var i e d, one can access various 
strongly interacting final states. On these figures, crosses + give the resuft of full numerical calculations, while the curves are 
our best fits using a Landau-Zener model. The functional forms for each fit are given in Table [n] (a) BBF L = + . (b) BBFF 
L — + . (c) BBFF L = 1~. Curves are fit from numerical result and "+" is numerical calculation. 



TABLE II: The final probability distribution of time evolution 
versus ramping speed is represented as a sequence of Landau- 
Zener transitions. Tij is defined as exp(— Xij/x)> where \ij 
is the fitted Landau-Zener parameter, and \ is the ramp ing 
speed H, A = The fitted XijS are listed in TabkjnTT 



BBF L = 


Pl (tnmci) 

p 2 (dimor + 1 atom) 

P3 ( 3 atoms) 


1 ~ T 12 

Tl 2 (l " T 23 ) 
T 12 T 23 


BBFF L — 


P1 (trimcr) 

P2 (dimer + 2 atoms) 
P3 (4 atoms ) 


(1 " T 12 )(l - T 13 ) 

Tl 2 (l " T 23 ) 

T 12 T 23 + C 1 — T 12) T 13 


BBFF L = 1 


p, (trimcr + 1 atom) 
p 2 (2 dimcrs) 
p 3 (1 diracr + 2 atoms) 
p 4 (4 atoms) 


(1 _ T 12 )(l - T 13 ) 

Tl 2 (l " T 23 )(l - T 24 ) 

(Tl 2 T 23 + (1 - T 12 )T 13 )(1 - T 34 ) 

T 12 T 23 T 34 + C 1 - T 'l2) T 13 T 34 + T 12( 1 - T 23) T 24 



speed. Moreover, when our few-body calculation are ap- 
plied to the many-body gas through a frequency rescal- 
ing, one might expect any any phase coherence at the 
few-body level could be washed out for the many body 
system, since there are multiple sources of decoherence 
in a many-body gas. 

C. P-matrix analysis 

The Landau-Zener parameters can alternatively be ex- 
tracted directly from the numerically calculated nonadi- 
abatic coupling P-matrix, as an alternative to the fit- 
ting procedure discussed above. Each P-matrix element 
connecting two low-lying states corresponds to a single 
important transition between state ^ and ^j, where 
the Landau-Zener parameter \ij can be extracted with- 
out ambiguity. The reasonable agreement with the fitted 
Landau-Zener model parameters suggest that those tran- 



sitions are adequately described by a two-level model, 
and that contributions from other levels' nonadiabatic 
coupling to other levels in the transition are negligi- 
ble. The dominant pathways are apparently isolated 
from others sufficiently far to justify this approximation. 
However, some P-matrix elements involving high-lying 
states have multiple-peaks, non-Lorentzian shapes for 
each peak, and apparent discontinuities in the explored 
region of ai lo /a sc . This phenomena could have several 
origins, the basic reason for its violation of Landau-Zener 
model validity is the non-isolated multi-channel nature of 
the transition, i.e. near the transition point, the nona- 
diabatic coupling between more than 2 levels at a time 
could be non-negligible. The following methods are used 
to address these issues. 

For the case of multiple-peaks in a give P-matrix cle- 
ment, the most important peak can be selected by consid- 
ering both the order and strength of the transition. For a 
non-Lorentzian shape, the most important peak is fitted 
as a Lorentzian form with two fitting parameters instead 
of using just one parameter as in the standard Landau- 
Zener model as reformulated by Clark [2Tj . The P-matrix 
elements for the most important transitions are shown 
in Fig. [3j along with the best fits and the full numeri- 
cal calculation. In a standard Landau-Zener model, the 
predictions extracted from the height and width of the 
P-matrix should coincide, but in our calculation not all 
the P-matrix elements have this feature because some ex- 
hibit non-Lorentzian shapes. For instance, consider Pi 3 
in the BBF system, for which the fitted Landau-Zener 
parameters extracted from our numerical time-evolution 
are seen to have values close to at least one of the two 
fitted P-matrix parameters or, they fall between the two 
alternative values fitted to either the height or the width 
of the assumed Lorentzian form of the Py . 

The trimer-dimer transition, corresponding to P12, is 
sharper than the atom-dimer transition, which suggests 
that adiabatic formation of trimer states requires a slower 
ramping speed than the formation of the dimer states. 
This fact is consistent with the result of our numerical 
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TABLE III: Comparison of Landau-Zener parameters ex- 
tracted directly from the nonadiabatic coupling P-Matrix cal- 
culation and from fitting the Landau-Zener model directly 
to our numerical time evolution. Here the terms height and 
width mean that the Landau-Zener parameters \ij below have 
been fitted to the height and width of the corresponding P- 
matrix elements. The incoherent transition functional form 
is listed in Table [TT| The number labeling of Xij is consistant 
with the labeling rule in Fig. [3] 
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time propagation. The Landau-Zener parameter charac- 
terizes quantitatively the typical ramping speed needed 
in order to stay in a trimer state or to ramp into either 
dimer states or atomic states. 
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3: Some important P-Matrix elements calculated as 
3X*j), in which X — 2&a. ( a ) BBF L = 0+. 
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FIG 

Pa = <*i 

(b) BBFF L = + .(c) BBFF L = 1 + . The inset of each figure 
contains a zoom in of the P-matrix elements except for P12. 
Both numerical calculated P-matrix elements and their fitted 
curve into Lorentzian form are shown. The number labeling 
of P-matrix (Pij) between different classes are: (a) 1:BBF, 2: 
BF+B, 3: B+B+F, (b) 1: BBF+F, 2: BF+B+F, 3: B+B+F, 

(c) 1: BBF+F, 2: BF+BF, 3: BF+B+F, 4: B+B+F+F. 



D. Molecule formation percentage 

Prediction of the molecule formation ratio as a func- 
tion of the inverse scattering length ramping speed is the 
strongest connection between our few-body calculation 
and the experimental many-body observations. Fig. [4] 
shows the molecule formation ratio in BBFF L = + and 
L = 1~ systems. Although the L = 1~ system has a 
higher upper limit for its molecule formation ratio, the 
formation rates can hardly be distinguished at fast ramp- 
ing speed. 

Various experiments have been carried out that ex- 
plored the maximum formation rate of hctcronuclear Fes- 
hbach molecules in recent years, including 6 Li 7 Li [2"6"ll27| . 
6 Li 23 Na [ 



40 K 87 Rb [H [291 Hi] and 6 Li 133 Cs [S]. 
The highest conversion through adiabatic magnetic field 
ramping across Fano-Feshbach resonance in 
given by JILA [HJ [T^] is around 36 percent of the mi- 
nority ( Rb). Our calculation gives the conversion rate 
as function of ramping speed in Figj4j By comparing the 
Bose-Fermi mixture experiment to former experiments, 
we notice that in Bose-Fermi mixtures the conversion 
rate never approach anywhere close to unity. The rea- 
son for this low conversion efficiency compared to the 
two-component fermion system and ultracold boson sys- 
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FIG. 4: The molecule formation percentage is depicted as 
a function of the dimensionless ramp speed for (a) L — + 
and (b) L — BBFF systems. Since the trimer and dimer 
were not distinguished experimentally, two separate analysis 
are presented here. In one analysis, a trimer is counted as a 
molecule (dashed) and in the other analysis the trimer formed 
is excluded from the molecules counted (solid) . 

tem [51 [HI [33] , which are described well by a stochastic 
phase-space paring model |32j , however this model is un- 
der debate. The predicted molecule formation rate in the 
stochastic phase-space pairing model is in general higher 
than the observed rate in magneto-association experi- 
ments [12] . In our few-body calculation, the possibility 
of trimer formation provides a way to potentially under- 
stand why the observed molecule formation rate might 
be reduced: the trimer is more deeply bound than the 
Feshbach molecules, and might decay rapidly into one 
deeply bound dimer if struck by another atom, or if it 
predissociates, in which case the "spin flip" method for 
detecting the dimer [H] would not be expected to detect 
a trimer bound state. 

Our few-body calculation could be directly applied to 
understand optical lattice experiments when the tunnel- 
ing between sites is small. Ideally, one boson and one 
fermion in one site is expected, and the presence of addi- 
tional bosons is likely to introduce three-body loss. How- 
ever, our calculation suggests that within an appropriate 
sweeping speed range, more than 2 particles on a single 
site could be directly converted into a molecular bound 
state. 



IV. DISCUSSION 

The Landau-Zener model yields a qualitative picture 
of the dynamical features of the trapped few-body sys- 
tem. It also offers some quantitative predictions for the 
transition probabilities and the characteristic range of 
ramp speeds that cause a change-over from a high rate of 
molecule formation to a low rate. However, the question 
of whether a Landau-Zener function is the appropriate 
functional form to describe the molecule formation frac- 
tion in a large system is still under debate [33, 34 . The 
Landau-Zener model (LZ) for three or four particles does 
not predict a single LZ function but rather a combination 
of different LZ terms that incoherently add up. Owing to 
its formulation in terms of the adiabatic eigenfunctions, 
the Landau-Zener model could identify each important 
transition into various possible final configurations rea- 



FIG. 5: Comparison of our numerical calculation carried out 
in the adiabatic representation of a two-level system (squares) 
with the numerical calculation in the diabatic representa- 
tion for a multilevel system (diamonds), and the analytical 
Landau-Zener model's prediction for the asymptotic behavior 
of the transition probability (solid line), (a) Probability to 
form a trimer state as a function of ramping speed a = 
in the two-level adiabatic calculation. The levels in included 
are the lowest BBF+F channel and lowest BF+B+F channel, 
(b) Cosine function of the phase accumulated during the tran- 
sition, where the phase accumulation is defined as the phase 
of the wave function minus the semiclassical phase. (This test 
is carried out in the BBFF L — case.) 



sonably well. Thus the "sequence of transitions" func- 
tional form appears to be consistent with our numerical 
simulations. But a two-level Landau-Zener model could 
not readily incorporate the effect of other nearby levels 
that may couple to one of the two levels, while the fi- 
nal probability is regarded as a sum of the probabilities 
for transitions into all levels in each class, which is a 
limitation of the Landau-Zener model for analyzing our 
numerical results. 

Several reasons might explain the discrepancies be- 
tween the P-matrix calculations and the fitted-time- 
evolution values for the Landau-Zener parameters. One 
reason for the discrepancies is that the Landau-Zener 
model as we have implemented it here does not include 
the coherent phase information of the time-evolution. 
Thus the phase coherence of the different indistinguish- 
able pathways could affect the final probability distri- 
bution, as in the Stueckelberg case. In the analytical 
solution of the two-level Landau-Zener model for arbi- 
trary initial conditions, the time evolution could be un- 
derstood as a "sudden" jump of phase plus the additional 
phase accumulated in the simple semiclassical time evo- 
lution asymptotically, since the semiclassical approxima- 
tion accurately describes the slowly varying phase. A 
numerical test shows that this "asymptotic" approxima- 
tion is only valid in the regime of very slow ramping speed 
in Fig. [5] For each transition, the prediction from the 
Landau-Zener model has better agreement with the ab- 
solute value of the transition probability than with the 
phase of the transition probability obtained in the nu- 
merical calculation. Evidently the LZ model, which pre- 
dicts the asymptotic behavior of transition probabilities 
between coupled levels, is more successful in predicting 
the probability information than the phase information. 

Another possible source of discrepancy is the fact that 
many levels have avoided crossings that are not suffi- 
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ciently well isolated, thus invalidating the Landau-Zcncr 
model. This error source differs from the first one in 
the sense that the second considers the effect of several 
crossings at the same time, while the first error source 
considers only isolated two-level crossings. 

In addition, considering our calculation's attractive 
Gaussian model potential, we point out that this po- 
tential corresponds to a specific three-body parameter 
which controls the trimer's size and binding energy. In 
our calculation, the three body parameter n is defined 
as \E trimer \ = ^ mer , where E tnmer is the trimer's 
binding energy at unitarity without trap, and fJ-trimer 
is the three body reduced mass of the trimer fitrimer = 
( y?' W 2 . For the 87 Rb and 40 K Bose-Fermi mix- 

V 7711+7712+7713 ; 

ture, fitrimer = 0.686, E tr imer = — 44.02/kj, where the 
dimensionless parameter has the value Ka^ Q = 5.49 in 
our units. For another boson-fermion interaction model 
potential, if it is within the same a^ Q /a sc range and has 
the same dimensionless three body parameter nciho, it 
will agree with the spectra and dynamics from our calcu- 
lation near unitarity (unless it is too far from unitarity, 
in which case the finite range correction of the potential 
is not negligible.) 

V. CONCLUSION 

Three and four-body problems remain fundamental 
and challenging. This study has presented an accu- 



rate numerical solution of the spectrum and dynamics of 
three and four-body Bose-Fermi mixtures in a trap in the 
vicinity of a Fano-Feshbach resonance. Even though the 
spectrum presents a rich structure of avoided crossings, 
a simple Landau-Zener model approximately describes 
the dynamics of unidirectional ramps and clarifies the 
relevant characteristic ramping speeds. The spectrum 
and dynamics are of immediate relevance to optical lat- 
tice experiments. These would allow access to physics 
that cannot be probed in a two-body system, such as 
trimer-atom interactions and dimer-dimer interactions. 
Also, the system of two bosons and two spin-polarized 
fcrmions in a trap exhibits many of the ingredients of 
the Fano-Feshbach resonances that form fermionic dipolc 
molecules. In that sense, our results provide a few-body 
perspective on many-body Bose-Fermi mixture experi- 
ments carried out at ultracold temperatures. 
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